/******************************************************************************/
/*                                                                            */
/*  Copyright (c) 2009: ngi-central.org                                       */
/*  All rights reserved.                                                      */
/*                                                                            */
/******************************************************************************/


/// \file
/// This file implements the algorithm for moments calculation.
//!
//! More information about calculation of moments can be found at 
//! \ref moments_algorithm.

#ifndef __sgl_moments_h__
#define __sgl_moments_h__

#include <sgl/global_settings.h>
#include <sgl/conditions.h>
#include <sgl/polygon.h>
#include <sgl/unnormalized_moments.h>

namespace sgl {


/******************************************************************************/
/*                                                                            */
/* moments                                                                    */
/*                                                                            */
/******************************************************************************/

template <
	/// The geometric primitive. 
	typename E>
/// Get the moments.
//!
//! This function calculates the moments of a geometric primitive.
//!
//! /return The moments of the geometric primitive.
inline sgl::unnormalized_moments<typename E::value_type> moments(
	/// The primitive.
	E const& primitive)
{
	sgl::unnormalized_moments<E::value_type> m;
	error_type_does_not_support_calculation_of_moments(); (*) // see below
	return m;
}

// (*) If you get an error in this line, this means that your type does not
//     support the calculation of moments. However, if you think it should 
//     support the calculation of moments, you should provide an overload of the
//     moments algorithm, which is able to work with your type.

template <
	/// The element type. 
	typename T>
/// Get the moments.
//!
//! This function calculates the moments of a box.
//!
//! /return The moments of the box.
inline sgl::unnormalized_moments<T> moments(
	/// The box.
	sgl::box<2,T> const& b)
{
	return sgl::moments(sgl::polygon<2,T>(b));
}

template <
	/// The element type. 
	typename T>
/// Get the moments.
//!
//! This function calculates the moments of a polygon.
//!
//! /return The moments of the polygon.
inline sgl::unnormalized_moments<T> moments(
	/// The polygon.
	sgl::polygon<2,T> const& poly)
{
	sgl::unnormalized_moments<T> m;

	for (size_t i=0; i!=poly.size(); ++i)
	{
		size_t i_minus_1 = 0==i ? poly.size()-1 : i-1;
		T common_term = poly[i_minus_1].x()*poly[i].y() - 
			poly[i].x()*poly[i_minus_1].y();

		m.m00() += common_term;

		m.m10() += common_term * (poly[i_minus_1].x() + poly[i].x());
		m.m01() += common_term * (poly[i_minus_1].y() + poly[i].y());

		m.m20() += common_term * (
			poly[i_minus_1].x()*poly[i_minus_1].x() + 
			poly[i_minus_1].x()*poly[i].x() +
			poly[i].x()*poly[i].x());
		m.m11() += common_term * (
			T(2)*poly[i_minus_1].x()*poly[i_minus_1].y() + 
			poly[i_minus_1].x()*poly[i].y() +
			poly[i].x()*poly[i_minus_1].y() +
			T(2)*poly[i].x()*poly[i].y());
		m.m02() += common_term * (
			poly[i_minus_1].y()*poly[i_minus_1].y() + 
			poly[i_minus_1].y()*poly[i].y() +
			poly[i].y()*poly[i].y());

		m.m30() += common_term * (
			poly[i_minus_1].x()*poly[i_minus_1].x()*poly[i_minus_1].x() + 
			poly[i_minus_1].x()*poly[i_minus_1].x()*poly[i].x() +
			poly[i_minus_1].x()*poly[i].x()*poly[i].x() +
			poly[i].x()*poly[i].x()*poly[i].x());
		m.m21() += common_term * (
			T(3)*poly[i_minus_1].x()*poly[i_minus_1].x()*poly[i_minus_1].y() + 
			poly[i_minus_1].x()*poly[i_minus_1].x()*poly[i].y() +
			T(2)*poly[i_minus_1].x()*poly[i].x()*poly[i_minus_1].y() + 
			T(2)*poly[i_minus_1].x()*poly[i].x()*poly[i].y() + 
			poly[i].x()*poly[i].x()*poly[i_minus_1].y() +
			T(3)*poly[i].x()*poly[i].x()*poly[i].y());
		m.m12() += common_term * (
			T(3)*poly[i_minus_1].y()*poly[i_minus_1].y()*poly[i_minus_1].x() + 
			poly[i_minus_1].y()*poly[i_minus_1].y()*poly[i].x() +
			T(2)*poly[i_minus_1].y()*poly[i].y()*poly[i_minus_1].x() + 
			T(2)*poly[i_minus_1].y()*poly[i].y()*poly[i].x() + 
			poly[i].y()*poly[i].y()*poly[i_minus_1].x() +
			T(3)*poly[i].y()*poly[i].y()*poly[i].x());
		m.m03() += common_term * (
			poly[i_minus_1].y()*poly[i_minus_1].y()*poly[i_minus_1].y() + 
			poly[i_minus_1].y()*poly[i_minus_1].y()*poly[i].y() +
			poly[i_minus_1].y()*poly[i].y()*poly[i].y() +
			poly[i].y()*poly[i].y()*poly[i].y());
	}

	m.m00() /= T(2);
	m.m10() /= T(6);
	m.m01() /= T(6);

	m.m20() /= T(12);
	m.m11() /= T(24);
	m.m02() /= T(12);

	m.m30() /= T(20);
	m.m12() /= T(36);
	m.m21() /= T(36);
	m.m03() /= T(20);

	return m;
}


} // namespace sgl

#endif // __sgl_moments_h__
